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ABSTRACT 


The finite element method is widely used to find the stress fields 
caused by external loads and temperature changes in composite materials. 
This report presents the first of a set of systematic studies of the rate 
of convergence and effect of singularities on such solutions. The problem 
considered is that of a finite crack in an infinite medium under anti- 
plane shear load. For this problem, it is shown that the nodal force at 
the trip of the crack accurately gives the order of the singularity, that 
energy release methods can give the strength to better than 1% with 
element size 1/10 the crack length and that nodal forces give a much 
better estimate of the stress field than do the elements themselves. 

This study is being extended to inplane shear, uniform tension and 
composite materials as well as finite bodies of rectangular form. 
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I. INTRODUCTION 

The finite element method is widely used to find stress fields in 
non-homogeneous elastic materials, such as composites. For fracture 
analysis the most important information is usually the strength and 
nature of the singularity at the tip of a crack or at an interface be- 
tween two laminae [e.g., 1,2]. Since the stress and strain fields are 
non-analytic at such singularities, it is difficult to improve the 
accuracy by the accepted procedure of averaging the stress or strain 
values in contiguous elements [3]. Pian and Tong [4] have shown that a 
minimum potential energy formulation gives nodal displacements which 
converge uniformly to the continuum solution for the displacement field. 
They and Perks [5,6] have discussed the convergence to the stress field 
in the vicinity of a singularity. The current method of choice is to 
fit the stress field near the tip by a region with a very fine mesh [7] 
or by using a special singular element [7,8]. The first method requires 
a large number of points and the second requires a prior knowledge of 
the nature of the singularity. 

From Pian's work, it would appear that more information can be 
had from the nodal forces than has so far been exploited for static 
solutions; some work has already been done on finding energy release 
rates for both dynamic [9,10] and static problems [11,14]. 

A finite crack in an infinite elastic solid loaded at infinity in 
tension (Type I), plane shear (Type II) or antiplane shear (Type III) 
are useful subjects for investigation of convergence because they are 
conveniently modeled with rectangular elements, the resulting matrix 
equations may be solved with difference equation techniques, and the 
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solutions are well known for the homogeneous case. This report gives 
extensive results for antiplane shear; work continues on the other two 
loadings and on laminates. 


II. NUMERICAL SOLUTION 

A finite element solution for strain field in an infinite solid 
with an infinite plane crack of width 2c subjected to uniform antiplane 
shear at infinity is a set of nodal displacements w- • with the indices 
i,j going from -« to +<». If we put i = j = 0 at the center of the crack, 

then w . . = -w. . and w. . = w- • so that we only need the solution 

”1 I ij 1 >“J ' 

for non-negative i and j. Solution details are given in Appendix I 

where it is shown that, for the residual problem, the nodal forces on 
the crack are f Q j = -1, -N < j < +N and the displacements w Q j = 0, 

1 N | ^ |j|. We intend to investigate the effect of mesh size on the con- 
vergence of the solution; to do so, it is convenient to convert the 
above problem into the finite equivalent of a boundary integral problem 
by finding the finite equivalent of the Green's Function (N = 1) and 
then using superposition. 

Pian shows [4] that the displacement field defined by the nodal 
displacements and linear or bilinear element displacement fields 
(Appendix I) gives a total potential energy greater than that in the 
continuum if the boundary nodes are loaded by forces equal to 


oj 


f j+l 

v 

j-1 


w(5hU)d? 


( 1 ) 


where w(?) is the (linearly varying) constrained displacement field. 
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For this to be valid, we must use the force boundary conditions of 
equation I -1 3b . In Appendix I, it is shown that the residual problem 
has a constant shear stress on the crack, so that the elements of f 0 are 
all -1. We get the solution in two steps; the first step is to find the 
displacement field for w Q0 = 1, w ■ = 0, j f 0. From equations 1-14 
and 1-17 in Appendix I: 


w 




N-l . 

c q /2 + z c jCOS ( j kn/N ) + ( -1 ) j c n /2 
J ~ 1 


( 2 ) 


If we now let N -> <*> 



f 1 

C(5) COSirj^ d£ 

J o 


(3) 


where k-rr/N -► n£ in equations 1-16 and 1 1 -4. From equation 1 1-2, it can 
be seen that 

■1 . 

c 1 ^) cosirj£ d 5 (3) 

0 

so that the nodal displacement field can be found by numerical integra- 
tion; we only need w^ for our work. f Qo is now found from equation 

1-6 or 1-7; if we divide each element of wt by f we have the first 

1 J oo 

row of the displacement field for a unit force at i=0,j=0 and f is 
found from equation I - 1 3b . The linear (triangular element) and bilinear 
(rectangular) element values for f^QQ a 9 ree< ^ within 4 x 10~ 10 (0.01%) 
and agreed with the continuum value for - 1 (Equation 1-3) within 
5 x 10'" 9 . 

For crack length of 2N elements, we may write: 

„ * (2N) . f (2N) 

I 0 0 



( 4 ) 
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where f^ 2 ^ is a column of l's, 2N-1 long and D-j is a (2N-1 ) x ( 2N- 1 ) 
matrix each row of which is f ' , found from equation I - 1 3b , with the 
f term on the main diagonal. Equation 4 can be reduced to N equations, 


since w 


oj 


w • , and solved for w . The nodal forces for N < j are now 

0 5 ~ J 0 


found by extending D-| into an( i writing 

D„ W < 2N > 


'2 a 


A 

= f 


where the elements of f are f . , j >_ N. 

o 0 > J 


III. RESULTS 

Using the bilinear rectangular elements, the nodal force at j = N+l 
is about 7.5% lower than t/t q at x = N/N+l . For a constant x, this error 
decreases almost linearly with 1/N so that linear extrapolation gives 
remarkably accurate results. For example, the finite element value for 
x = 1.25 c is 1.5962 for the 8 element solution and 1.6302 for 16. Linear 
extrapolation to 1/M = 0 gives 1.6643; the continuum value is 1.6667, an 
error of 0.14%. Figure 1A shows T y Z / T 0 vs r/c as given by equations 
1-3 and 1-4. It would be difficult to distinguish the finite element 
values from the exact curve for j > N+l but F^ is shown for N from 2 to 
64. The trend for the constant strain (linear) elements is about the same 
but with about twice the error. 

The above comparison shows, as is well known, that the nodal dis- 
placement values converge to the continuous displacement function at a 
rate appropriate to the order of the equivalent finite difference 
algorithm and that the stresses at any analytic point may be found by 
numerical differentiation of the displacement field to the same order 
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by using Hooke's law. Averaging the "stresses" in contiguous constant 
strain elements gives less accurate values unless the weighting points 
are carefully selected, e.g., the Gaussian integration points [3]. 

One would expect [4] the singular force at j = N to be given by the 
weighted integral of the stresses on the y = 0 face from j = N-l to N+l . 
If we consider the original (not the residual) problem, x = 0, j < N. 

Then using equation 1-3 and a linear weighting function 


F 


N 


1+h 


x(l - x/h) 


dx 



(1 + o(h)) 


= 0.9428 Zh 


where h = 1/N. The continuum residual problem converted to finite element 
form differs from the residual finite element problem by adding 0.5 h to 

F r 

Since the displacement field solution is not analytic at the 
singularity, equations 1-9 and I -1 1 don't apply. However, we are inte- 
grating only on a line where w = o and ^ = 0 so that we may let w Q j = 0 

9X 

in these equations; if we use only w^ w-| ^ and w-j ^ we should get 
, o(h 4 ). These and several other less accurate (o(h^)) approximations 
were used and they all converged very rapidly to c^/h. The points for 
the bilinear case are shown in Figure lb. Rather than attempting to find 
the slope from the log-log plot, it is better to extrapolate F^/Zh vs h 
to h = 0. These curves are shown in figure 2 where we compare the nodal 
force at i = 0, j = N using equation 1-11 modified by letting w . e 0, 

0 5 J 

the linear approximation for I— and the second order approximation in which 
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we use 




These algorithms are: 


F N " 3 * W 1,N-1 + W 1,N + ”l,N+P * Eq 1-11 mo<i * 

F N = W 1,N = If + °< h2 > <« 

f n “ - 1 (w i,n-i + 4w i,n + w Tn+i> * a7 + 0(h4) (7) 


All of the curves converge to a square root singularity but the bilinear 
ones do so in a linear fashion which makes extrapolation accurate for 
crack lengths as coarse as 16 elements (N = 8). Unfortunately, they don't 
extrapolate to the correct coefficient, which shows that higher order 
terms must be included even in the limit. Table I shows the rate of con- 
vergence for to its value at h = 0 using equation 7 for both the linear 
and bilinear case and to the square root singularity for equations 6 and 7. 

The external work done is another measure of the strength of the 
singularity. The external work, equal to the strain energy, is 


tt/ 2 = U(c = 1) 


,1 

w(x)i yz (x)dx 


* 


N- 

j=4n-i) 


w 0J .(x) or Simpson's 


( 8 ) 


because F . = 1 . As shown in Table 1, this converges as rapidly as F^, and 
converges to the correct value. The columns headed 6 are the result of 
linear extrapolation using the entry and the previous entry; thus the 
extrapolation of U from the bilinear solution using 2 and 4 elements for 
the half crack length gives it/4 with an error of 0.3%. One may also 
approximate the energy release rate from 2 • (w^ • F^/2). This should 
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converge to 



(9) 


The % error for the linear and bilinear case and for the extrapolated values 
are shown in the last four columns of Table I. These results are better 
than Rybicki [14] indicates. 

Since the values for w Q agreed within 5 x 10 in the linear and 

bilinear cases and the small table top computer used in this work took 

about one minute per integration (equation 3), the values for j > 50 from 

the linear case were used for the bilinear solutions. This small in- 

-4 

consistency accounts for the non-uniform convergence at levels of 10 
for N = 64 in Table I. 


IV. CONCLUSIONS 

It has been shown that the singularity at a crack tip under 
anti plane shear (Type III) can be accurately modeled with a relatively 
coarse mesh near the tip if one uses the energy release from crack exten- 
sion. The results are remarkably accurate if one gets solutions for two 
or three element sizes and extrapolates to zero size (Table I). 

Curve fitting to the stress field requires a very fine mesh both 
because the slope of the continuum solution departs rapidly from that 
of the leading term and because the finite element solution is oscillating 
at the first two or three nodal points (see also Wilson [7]). 

The nodal force at the crack tip as a function of mesh size is a 


very accurate measure of the order of the singularity but may converge 
to the wrong strength if higher order terms are not included. Even 
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where the displacement field is analytic, the nodal forces give a better 
approximation to the stress field than does element averaging. 

The methods used in this report are being extended to Types I and 
II loading and laminated material and could be used for finite bodies 
of rectangular shape. They are economical with computer time since by 
using boundary element techniques, only small matrices need be inverted 
to get finite element solutions for large (or infinite) arrays. The 
results presented here were obtained on a small (64K bytes) table top 
computer. 
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APPENDIX I - FINITE ELEMENT PROBLEM FORMULATION 


Our problem is that of a crack of length -c £ x £ c in the y = 0 

plane, extending to ± °° in the z direction of the infinite solid with a 

2 2 

uniform shear stress x^ z = x Q at x + y -* ® (anti plane, Type III). This 
is the same problem as the edge crack of length c in the semi-infinite 
solid, x >_ 0. The continuum solution is given by Bentham and Koiter 
[12, Sect 3.3]. For antiplane shear: 


T xz '■ 1 T yz = 2<t> ' U) 1 


t = x + iy 


Gw = <|>U) + 4>(s) 
i = /T 


i-l 


For our geometry: 


♦to - - i i x o - c 2 ) 1 ' 2 


Along the x axis 


T = T 

yz o 


2 2 
x - c 


1/2 


x 2 >c 2 


2 2 

There is a singularity of order 1/2 at x = c . For > 0, let 
x = c + r,ct = £- 

Then x = T (2a) ^ 2 (1 + i- a + o(a 2 )) 

yz o 4 ' 

and K ttt = -r c^ 2 

1 1 1 0 


1-2 


1-3 


1-4 


The finite element problem is more easily handled if we subtract a 
uniform stress field from the above, leaving a boundary condition of: 
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T = - T 

yz o 


-> 0 


y = o 

2 2 
x + y ■> « 


2 2 
x < c 


1-5 


We use a square grid of nodal points at x = j, y = i subjected to 
a shear stress x = 1 on the negative face i = 0, -c < j < c and let G = 1 
By symmetry: 


W -i >j " W i,j 


W i,-j " W i > j 


The basic element is the square of side 1. By minimizing the strain 
energy [3] we can show that the nodal force at point i , j of the element(s) 
bounded by i, i+1 , j and j+1 is: 


F. • = w. . 
iJ 


2 w i+l,j ' 2 w i,j+l 


1-6 


for a uniform strain field in the two triangular elements bounded either 
by the line (i,j)(i+l ,j+l) or (i+1 ,j)(i,j+l). It is interesting and 
useful that for antiplane shear the direction of the diagonal does not 
affect equation 1-6. For the bilinear displacement field: 

w = w^j (l-c)O-n) + w i 4 -i 5 j cO _ n) 

+ w iJ+1 0-c)n + w i+1J+1 cn 
the corresponding equation is: 

F i j = I W i j " \ w i+l , j “ 1 w i , j+1 " 1 w i+l ,j+l 1-7 

These elements are now assembled to give a zero nodal force for i + 0. 
Using equation 1-6 the field equations are 

F ij ” 0 = - w i-l.j - + 4w i,o ' w i.jH - "1+1 .j 


1-8 



n 


and the nodal force on the negative y face is 


F i7 = + 2w i,j - - w i+i , j 


1-9 


Using equation 1-7, they are 


and 


0 = ' * Vl,j + w i-l,j + l> 


1 


1 


+ <■ i + 1 w ij - i w i , j+i > 

‘ I (w i+l,j-l + w i+l,j + w i+l,j+l ) 


1-10 


F ij = “ 3 (w i-l,j-l + W i-l,j + W i-l,j + l } 

+ ^ ( "i w i,j-i + f w i,j ~ i w i,j+i } 


i-n 


Equation 1-8 is the finite difference equation for the Laplacian 

O 

good to o(h ); using the governing partial differential equation 


2 2 
8 W 9 W 


3X 


9y 


1-12 


2 

we find that equation 1-9 gives the first derivative, w , o(h ), while 

w i t " w i 4 .i -i 1S o ( h ) . Normally we can use i (w. •, • - w.,, •) which is 
o 

o(h ) but not at a stressed surface. This same result does not hold for 

/ 4\ 

equations I -10 and 11 but we can show that the system converges o(h ). 
The corresponding matrix equation is 


- Bw.j -j + AvT - BWj + -j - 0 

I - 1 3a 

E??0 - * \ Aw 0 - Bw, = f 0 

I - 1 3b 


where w!j is the column vector of the displacements of the nodes in i* 
row. Using the notation of Appendix II, the values of a and b in equa- 
tion II-l are those in Table I. 
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Equation 

1-8 

1-9 

1-10 

i-n 

Matrix 

A 

B 

E F 

A 

B 

E 

F 

a 

4 

-1 

2 -1 

8/3 

- 1/3 

4/3 

- 1/3 

b 

-1 

0 

- 1/2 0 

- 1/3 

- 1/3 

- 1/6 

- 1/3 


Table 1-1 


We may solve the matrix difference equation I -13a, by expanding 
w. in the eigenvectors of Appendix II. Using equation I I -7, if we 
substitute: 


-rr +* \ 

w.j = X q. 

A = X A a Y • 

B = X A b I , 


1-14 


into equation I - 1 3a and use the last of equations II-2, the equations for 
'q. are decoupled: 


- a b Vl + a a q i - A B q i + 1 = 0 


1-15 


If we now let: 

then 

and 

where 


q ik “ c k q i-l ,k 
X Bk c k + A Ak " X Bk c k ~ 0 

c k = B k J ■V TT 


6 k x Ak /21 Bk 


1-16 


and Agj^ are given by equation 1 1 -4 and a and b are shown in table 1-1. 

In order to satisfy the condition at infinity, we take the + or - sign 

? -1 
to make c < 1. c^ ■> 1/2 B as B ->-0. 
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If w Q0 = 1, w Q j = 0, j f 0, from equations I - 1 3a and 1-16 

A Ak q ik A Bk C| 2k A Bk x ok w oo A Bk w oo 
^ A Ak ” c A Bk^ q ik = A Bk w oo 


’ik = C k w oo 
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APPENDIX II - FACTORING TRIDIAGONAL MATRICES 

The matrices we are concerned with all have one constant, a, on 
the main diagonal and another, b, on the first sub and superdiagonal 
with all other elements zero. By symmetry, they extend from -N to +N 
and are (2N+1) x (2N+1) square. In order to get symmetrical reflection 
at ± N, we make the second element of the first row and the next to the 
last element of the last row equal to 2b which gives the following form. 

a 2b 0 0 . . 

b a b 0 . 

0 b a b . 

"6 b *a”b’*6” 

. 0 2b a 

The associated eigenvalue problem is: 

RX = X A 
YR = AY 

f 1 1-2 

normalized to Y X = I 

so that R = XAY and YRX. = A 

where A is the diagonal matrix of the eigenvalues, X is the square matrix 
of the eigenvector columns and Y is the adjoint square matrix of the 
eigenvector rows. The difference equations are: 

b x i-l,k + a x i,k + b x i+l,k = X k x i,k 
a X -N,k + 2b X -N+1 ,k = A k X -N,k 
2b X N-1 ,k + a x N,k = X k x N,k 
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The matrix has 2N+1 roots: 

A k = a + 2b cos (kir/2N) k = 0,1,...,2N+1 II-4 

The corresponding column eigenvectors are: 

X ik = cos (ik /2N) i = -N, . . . ,0, . . . ,N 1 1 -5 

The adjoint row eigenvectors are: 

Y ki ■ <’ -r'l.tNX 1 -K N> k)*ik /N "- 6 

where <5.. is the Kronecker 6. 

’ J 

The antisymmetric tridiagonal matrix found by replacing the 2b 
terms in equation I I -1 with 0, has the same eigenvalues, but sines for 
the eigenvectors (equation I I -5 ) instead of cosines. 
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